Setup

library(repro)
# load packages from yaml header
automate_load_packages()
# include external scripts
automate_load_scripts()

# load data
intern  <- automate_load_data(intern, read.csv, stringsAsFactors = T)
medic    <- automate_load_data(medic, read.csv, stringsAsFactors = T)
PAV      <- automate_load_data(PAV, read.csv, stringsAsFactors = T)
INST     <- automate_load_data(INST, read.csv, stringsAsFactors = T)
PIT      <- automate_load_data(PIT, read.csv, stringsAsFactors = T)
HED      <- automate_load_data(HED, read.csv, stringsAsFactors = T)
HED_fMRI <- automate_load_data(HED_fMRI, read.csv, stringsAsFactors = T)




## we recommend running this is a fresh R session or restarting your current session
#install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
#install_cmdstan()


# check_git(); check_make(); check_docker() #check if installed


sessio = session_info(); #opts_chunk$set(echo = F, message=F, warning=F) # set echo F for all


#May I suggest running `repro::automate()`? 

#This will create a `Dockerfile` & `Makefile` based on every RMarkdown in this folder and the special yamls in them. date: "`r format(Sys.time(), '%d %B, %Y')`" 

#add ENV DEBIAN_FRONTEND=noninteractive to DOCKERFILE

This file was automatically created via the Repro package (version 0.1.0) using R version 4.0.4 (2021-02-15)

niter = 5000; warm = 1000; chains = 4; cores = 4; nsim = 10000 # number of iterations (to change if you want to quick check and warmups (BUT chains and BF might be really unstable if you have less than 20'000 iter (4x5000) ) #or also parallel::detectCores()/2)
options(scipen = 666, warn=-1, contrasts=c("contr.sum","contr.poly"), mc.cores = cores)  #remove scientific notation # remove warnings #set contrasts to sum ! #remove scientific notation # remove warnings #set contrasts to sum !
 #cl = parallel::detectCores()/2
set.seed(666) #set random seed
control = lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')) #set "better" lmer optimizer #nolimit # yoloptimizer
#emm_options(pbkrtest.limit = 8000) #increase repetitions limit for frequentist stats

source('R/plots.R', echo=F)# plot specification
source('R/utils.R', echo=F)# useful functions

panderOptions('knitr.auto.asis', FALSE) #remove auto styling

labels <- c("-1" = "Pre", "1" = "Post")#for plots

# Look at R/clean.R (listed in the YAML) which does all the preprocessing for more info


# If you are unsure weather or not you have `git` `make` & `docker`.
# check_git()
# check_make()
# check_docker()

Description

Stan uses Hamiltonian Monte Carlo (HMC) to explore the target distribution — the posterior defined by a Stan program + data — by simulating the evolution of a Hamiltonian system.

Demographics

egltable(c("BMI", "AGE", "GENDER"), 
  g = "INTERVENTION", data = df, strict = FALSE) %>%
  kbl(caption ="Summary statistics", digits = 2) %>%
  kable_styling(latex_options = "HOLD_position", position = "center", full_width = F) %>%
  row_spec(0,bold=T,align='c')
Summary statistics
Placebo M (SD)/N (%) Liraglutide M (SD)/N (%) Test
BMI 35.08 (2.96) 35.74 (2.93) t(df=47) = -0.78, p = .438, d = 0.22
AGE 40.27 (13.74) 38.61 (11.72) t(df=47) = 0.45, p = .653, d = 0.13
GENDER Chi-square = 0.35, df = 1, p = .556, Phi = 0.08
Men 10 (38.5) 7 (30.4)
Women 16 (61.5) 16 (69.6)

Biomedical data

Variable Selection

# Boxplot of biomedical variables per group
ggplot(med)+
  geom_boxplot(aes(INTERVENTION, n))+
  facet_wrap(~feature, scales = "free")+
  labs(title = "")+
   theme_minimal()+
  theme(axis.title = element_text()) + 
  ylab("Biomedical predictor's value (scaled)") + 
  xlab('')
Box-plot of all biomedical predictors per intervention.

Box-plot of all biomedical predictors per intervention. 

# Plot correlogram of numeric variables
#pairs(~., data = df[,8:19], main = "Scatterplot Matrix of variables")
#corrplot(cor(df[,8:19], use="pairwise.complete.obs"), type="lower")


Recursive Feature Eliminations

#1) with Recursive Feature Eliminations (CARET)

sizes = 1:length(dfmed[c(-1)]); len = length(sizes)
seeds <- vector(mode = "list", length = len)
for(i in 1:(len-1)) seeds[[i]]<- sample.int(n=nsim, size = length(sizes)+1)
# for the last model
seeds[[len]]<-sample.int(nsim, 1)

RFEcontrol <- rfeControl(functions=rfFuncs, method="cv", number=10, seeds= seeds) # control options

rfeResults = rfe(x = dfmed[c(-1)], y = dfmed$intervention, sizes=sizes, rfeControl=RFEcontrol)
predictors(rfeResults)
## [1] "BMI_diff"    "BW_diff"     "reelin_diff" "GLP_diff"
plot(rfeResults, type=c("g", "o")) # look for the "elbow"

#if we agree that BMI, Body Weight and Waist Circumference actually measure the same thing, there only 4 other variables that are "useful" to separate the two groups :
#Reelin and GLP

Mediation analysis

#parallel:::setDefaultClusterOptions(setup_strategy = "sequential")

med_res <- intmed::mediate(y = "weightLoss", med = c("GLP_diff" ,  "reelin_diff"),  treat = "intervention", ymodel = "regression", mmodel = c("regression", "regression"), treat_lv = 1, control_lv = 0, incint = TRUE, inc_mmint = FALSE, conf.level = 0.95, data = df , sim = nsim, complete_analysis = TRUE, digits = 3,  summary_report=F) #c = c("age","gender"),

table <- list.clean(readHTMLTable("res.html"), fun = is.null, recursive = FALSE); table = table[3]$`NULL`[1:6,]
#pander(table)
colnames(table)[4] = "p"; table$p = as.numeric(table$p); table$p = ifelse(as.numeric(table$p) < 0.05,paste("<span style=\" font-weight: bold; \" >" ,sprintf("%.3f",table$p), "</span>"),  paste("<span>" ,sprintf("%.3f",table$p), "</span>")) # highlight p < 0.05
table$p = ifelse(table$p == '<span style=" font-weight: bold; " > 0.000 </span>', "<span style=\" font-weight: bold;    \" >\u003C 0.001</span>", table$p)


table[-c(3:4),] %>%
  kbl(caption ="Mediation Analysis: DV = Weight loss, IV = Intervention", escape=F) %>%
  kable_styling(latex_options = "HOLD_position", position = "center", full_width = F)
Mediation Analysis: DV = Weight loss, IV = Intervention
Estimates 95% CI p
1 Indirect effect mediated through GLP_diff 0.017 (-0.278, 0.327) 0.914
2 Indirect effect mediated through reelin_diff 0.048 (-0.105, 0.248) 0.543
5 Direct effect 1.425*** (0.968, 1.884) < 0.001
6 Total effect 1.489*** (1.024, 1.981) < 0.001 
#just for plot purpose (because not the same exact analysis)
medi =  psych::mediate(weightLoss ~ intervention + (GLP_diff) + (reelin_diff), data = df, n.iter = nsim, plot=F); psych::mediate.diagram(medi, show.c=FALSE)

Weight Loss

# Model
mf = formula(weightLoss ~ intervention + gender + age + GLP_diff + reelin_diff  + (1|id))

# -------------------------------------- Bayesian Regression Models using Stan -------------------------------------
# STAN is a probabilistic programming language that allows you to get full Bayesian statistical inference with MCMC sampling.


bmod_full = brm(mf, data=df, family = gaussian, prior = c(prior(normal(0, 3), class = "b", coef = ""), prior(normal(0, 100), class = "Intercept", coef = "")), sample_prior = T, save_pars = save_pars(all = TRUE), chains = chains,  iter = niter, warmup = warm, seed = 123, inits = 1, backend = 'cmdstanr', threads = threading(4), control = list(adapt_delta = 0.99)) # a lot to unwind here..  1) Generic informative prior around 0 for fixed effects and weak prior for the intercept 2) we need to sample priors and save parameters for computing BF 3)larger step size  
#lmer to compare
fmod_full = lm(update(mf, ~.- (1|id)) , data = df)
## plot population-level effects posterior distributions and chain sampling

param = mcmc_plot(object = bmod_full, pars =c("b_.*en", "b_.*ge"), type ="areas")

trace = mcmc_plot(object = bmod_full, pars =c("b_.*en", "b_.*ge"), type ="trace")

#check assumptions
var_group = pp_check(bmod_full, type = "stat_grouped", group = "intervention", binwidth = 0.1, nsamples = NULL) #equality of variance between groups

rep_fit = pp_check(bmod_full, nsamples = 100) # check response fit

error = pp_check(bmod_full, type ="error_scatter_avg", nsamples = NULL) # check good alignment between model and data, and no obvious pattern to the types of errors we are getting.


#Normality of errors
residuals <- residuals(bmod_full)[, 1]; res <- qqnorm(residuals, pch = 1, plot.it=F)

lmx = plot_model(fmod_full, type = "diag"); #diagnostic plots for lmer

#plot all
diagWEIGHT <- ggarrange(param, var_group, rep_fit, error, ncol = 2, nrow = 2)

annotate_figure(diagWEIGHT, top = text_grob("Diagnostic Plots", face = "bold", size = 10))

full_tab = describe_posterior(bmod_full,
                   estimate = "median", dispersion = T,
                   ci = .9, ci_method = "hdi",
                   bf_prior = bmod_full, diagnostic = "Rhat",  
                   test = c("p_direction", "bf"))
full_tab = filter(full_tab, Parameter %notin% c("b_Intercept", "prior_sigma"))



# -------------------------------------- FREQUENTIST STATS: AOV -----------------------------------------------

# formula = 'weightLoss ~ intervention + gender + age  + Error(id)'
# 
# 
# mdl.weight = aov_car(as.formula(formula), data = df, factorize = FALSE,  anova_table = list(correction = "GG", es = "pes"),type = "II")
# res = nice(mdl.weight, MSE = F); ref_grid(mdl.weight); res
# 
# PES.weight = pes_ci(as.formula(formula), data = df, conf.level = .90, factorize = FALSE, anova.type = "II", epsilon="none") ; 
# mdl.weight.emms = emmeans(mdl.weight, pairwise ~ intervention)


# res$p.value = as.numeric(res$p.value)
# res$p.value = ifelse(res$p.value < 0.05,paste("<span style=\" font-weight: bold; \" >" ,sprintf("%.3f",res$p.value), "</span>"),  paste("<span>" ,sprintf("%.3f",res$p.value), "</span>"))
# 
# res$F = unlist(str_split(gsub("[^0-9.,-]", "", res$F), ","));res$pes = unlist(str_split(gsub("[^0-9.,-]", "", res$pes), ","));
# res$`90% CI` = paste(sprintf("%.3f",PES.weight[,2]), "-", sprintf("%.3f",PES.weight[,3]))
# 
# res$p.value[1]= "<span style=\" font-weight: bold;    \" >\u003C 0.001</span>"
# colnames(res)[3:5] = c( paste("F(", res$df[1], ")", sep=""),"&eta;<sub>p</sub><sup>2</sup>", "p")
# res[c(1,4,6,3,5)]  %>% kbl(digits = 2, escape = F) %>%
#   kable_styling(latex_options = "striped", position = "center", full_width = F) 

# -------------------------------------- Regression table summary --------------------------------------

tab_model(bmod_full, show.p = F, show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 3, dv.labels = "Weight Loss (\u0394 BMI)", file = "tmp/temp1.html", transform = NULL,  rm.terms = "beta")

report = capture.output(sexit(bmod_full, ci=.9))
print(paste("Bayesian linear mixed model (estimated using MCMC sampling with " ,chains ," chains of", niter, " iterations and a warmup of", warm, ") to predict weightLoss with intervention, gender, age, GLP_diff and reelin_diff (formula: weightLoss ~ intervention + gender + age + GLP_diff + reelin_diff). The model included id as random effect (formula: ~1 | id). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 2.50) distributions"))
## [1] "Bayesian linear mixed model (estimated using MCMC sampling with  4  chains of 5000  iterations and a warmup of 1000 ) to predict weightLoss with intervention, gender, age, GLP_diff and reelin_diff (formula: weightLoss ~ intervention + gender + age + GLP_diff + reelin_diff). The model included id as random effect (formula: ~1 | id). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 2.50) distributions"
full_tab
## Summary of Posterior Distribution
## 
## Parameter    | Median |  MAD |        90% CI |     pd |  Rhat |     BF
## ----------------------------------------------------------------------
## intervention |   1.36 | 0.22 | [ 0.98, 1.71] |   100% | 1.000 | > 1000
## gender       |   0.16 | 0.19 | [-0.19, 0.48] | 78.53% | 1.001 |  0.090
## age          |  -0.29 | 0.20 | [-0.62, 0.01] | 93.47% | 1.002 |  0.205
## GLP_diff     |  -0.19 | 0.28 | [-0.65, 0.31] | 74.41% | 1.001 |  0.121
## reelin_diff  |  -0.04 | 0.29 | [-0.53, 0.42] | 55.42% | 1.000 |  0.100
report[4]
## [1] "- b_intervention (Median = 1.36, 90% CI [0.98, 1.71]) has a 100.00% probability of being positive (> 0), 100.00% of being significant (> 0.09), and 99.98% of being large (> 0.54)"
tables <- list.clean(readHTMLTable("tmp/temp1.html"), fun = is.null, recursive = FALSE)
tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
tables2 <- as.matrix(tables2) %>% as_tibble();


tables2[is.na(tables2)] <- "";  names(tables2) <- NULL; tmp = t(tables2[(length(full_tab$Parameter)+1):(length(full_tab$Parameter)+4),1:2]);
tmp[,4][1] = "R2"; tmp[,4][2] = gsub(".*/","",tmp[,4][2]); colnames(tmp) =  tmp[1,]
pander:: pander(tmp[2,])
## 
## ------------------------------------
##  ICC    N id   Observations    R2   
## ------ ------ -------------- -------
##  0.16    45         45        0.730 
## ------------------------------------

Plot Weight Loss

dfdraws = bmod_full %>%
    spread_draws(`b_intervention` )


HDI_weight = plotHDI( dfdraws$`b_intervention` , credMass = .90, binSize = 100, Title = "") + theme_bw() + html_theme  


dfdraws2 =  bmod_full %>%
     emmeans(~ intervention) %>%
     gather_emmeans_draws() 


pp = dfdraws2 %>%
    ggplot(aes(x = as.factor(intervention) , y = .value,  fill = as.factor(intervention))) +
    geom_abline(slope= 0, intercept=0, linetype = "dashed", color = "black") +
    stat_slab(.width = c(0.50, 0.9), position="dodge", alpha=0.5) +
    stat_pointinterval(.width = c(0.50, 0.9),position="dodge") +
    labs(x = "",y = "Weight Loss (\u0394 BMI)", title = "") + 
    scale_fill_manual(values=c("-1" = pal[1],"1"=pal[6]), guide="none") +
    scale_color_manual( values=c("-1" = pal[1],"1"=pal[6]), guide="none") +
    scale_x_discrete(labels=c("Placebo", "Liraglutide")) +
  #scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(400,550, by = 50)), limits = c(380,575)) +
    theme_bw()

plt_bayes_html = pp + html_theme 
plt_bayes = pp + averaged_theme 


figureWL <- ggarrange(plt_bayes_html, HDI_weight,
                    labels = c("A", "B"),
                    ncol = 2)
#figureWL




dfR <- summarySE(df, measurevar = "weightLoss",
                 groupvars = "INTERVENTION")

dfR$cond <- ifelse(dfR$INTERVENTION == "Liraglutide", -0.25, 0.25)
df$cond <- ifelse(df$INTERVENTION == "Liraglutide", -0.25, 0.25)
df <- df %>% mutate(condjit = jitter(as.numeric(cond), 0.3),
                    grouping = interaction(id, cond))

pp <- ggplot(df, aes(x = cond, y = weightLoss, 
                     fill = INTERVENTION, color = INTERVENTION)) +
  geom_hline(yintercept=0, linetype="dashed", size=0.4, alpha=0.8) +
  geom_flat_violin(scale = "count", trim = FALSE, alpha = .2, aes(fill = INTERVENTION, color = NA))+
  geom_point(aes(x = condjit), alpha = .3,) +
  geom_crossbar(data = dfR, aes(y = weightLoss, ymin=weightLoss-se, ymax=weightLoss+se), width = 0.2 , alpha = 0.1)+
  ylab('Weight loss (\u0394 BMI)')+xlab('')+
  scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(-2.5,7.5, by = 2.5)), limits = c(-3,8)) +
  scale_x_continuous(labels=c("Liraglutide", "Placebo"),breaks = c(-.25,.25), limits = c(-.5,.5)) +
  scale_fill_manual(values=c("Liraglutide"= pal[6], "Placebo"=  pal[1]), guide = 'none') +
  scale_color_manual(values=c("Liraglutide"= pal[6], "Placebo"=  pal[1]), guide = 'none') +
  theme_bw() 

plt = pp + averaged_theme 
pp + html_theme
Weight Loss by intervention.

Weight Loss by intervention. 



cairo_pdf('figures/Figure_WL_Lira.pdf')
print(plt)
dev.off()

cairo_pdf('figures/Figure_WL_Lira_Bayes.pdf')
print(figureWL)
dev.off()

Pavlvovian Conditioning Task

Latency

Latency = time to detect the target (ms) & condition = CS+ or CS-

# Model
mf = formula(RT ~ condition*intervention*session + age + gender + BMI_V1  + hungry + GLP_diff + reelin_diff + (condition*session|id) + (1|trialxcondition))


# -------------------------------------- Bayesian Regression Models using Stan -------------------------------------
# STAN is a probabilistic programming language that allows you to get full Bayesian statistical inference with MCMC sampling.
bmod_full = brm(mf, data=PAV, family = exgaussian, prior = c(prior(normal(0, 3), class = "b", coef = ""), prior(normal(0, 100), class = "Intercept", coef = "")), sample_prior = T, save_pars = save_pars(all = TRUE), chains = chains,  iter = niter, warmup = warm, seed = 123, inits = 1, backend = 'cmdstanr', threads = threading(4), max_treedepth=15) # a lot to unwind here..  1) Generic informative prior around 0 for fixed effects and weak prior for the intercept 2) we need to sample priors and save parameters for computing BF  3) account for the skewness of reaction times #family = "exgaussian"
#lmer to compare
fmod_full = lmer(mf , data = PAV, REML=F, control = control )
## plot population-level effects posterior distributions and chain sampling

param = mcmc_plot(object = bmod_full, pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="areas")

trace = mcmc_plot(object = bmod_full, pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="trace")


#check assumptions
var_group = pp_check(bmod_full, type = "stat_grouped", group = "intervention", binwidth = 1, nsamples = NULL) #equality of variance between groups

rep_fit = pp_check(bmod_full, nsamples = 10) # check response fit -> accounting for skewness, nice!

error = pp_check(bmod_full, type ="error_scatter_avg", nsamples = NULL) # check good alignment between model and data, and no obvious pattern to the types of errors we are getting.


#Normality of errors
residuals <-residuals(bmod_full)[, 1]; res <- qqnorm(residuals, pch = 1, plot.it=F)

lmx = plot_model(fmod_full, type = "diag"); #diagnostic plots for lmer

#plot all
diagRT <- ggarrange(param, var_group, rep_fit, error, ncol = 2, nrow = 2)

annotate_figure(diagRT, top = text_grob("Diagnostic Plots", face = "bold", size = 10))

full_tab = describe_posterior(bmod_full,
                   estimate = "median", dispersion = T,
                   ci = .9, ci_method = "hdi",
                   bf_prior = bmod_full, diagnostic = "Rhat",  
                   test = c("p_direction", "bf"))
full_tab = filter(full_tab, Parameter %notin% c("b_Intercept", "prior_beta", "prior_sigma"))



# -------------------------------------- FREQUENTIST STATS: LRT -----------------------------------------------
#measurevar = "RT"
# formula = "condition*intervention*session"
# COVA = "+ condition*age + gender + BMI_V1  + thirsty  + hungry +"
# random = " (condition*session|id) + (1|trialxcondition)"
# 
# #method LRT 
# model = mixed(as.formula(paste(measurevar, paste (formula, COVA, random), sep="~")), data = PAV, method = "LRT", control = control, REML = FALSE); model
# 
# table = nice(model);  ref_grid(model)  #triple check everything is centered at 0
# 
# 
# mod =  lmer(as.formula(paste(measurevar, paste (formula, COVA, random), sep="~")), data = PAV, control = control)

# tab_model(mod, show.p = T,show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 2, dv.labels = "Latency", emph.p = TRUE, file = "tmp/temp1.html")

# tables <- list.clean(readHTMLTable("tmp/temp1.html"), fun = is.null, recursive = FALSE)
# tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
# 
# tables2 <- as.matrix(tables2) %>% as_tibble()
# tables2[is.na(tables2)] <- ""
# tables3 = tables2[1:length(table$Effect),1:4] 
# tables3[5] = str_split(gsub("[^0-9.,-]", "", table[3]), ",")[[1]]; tables3[6] = as.numeric(str_split(gsub("[^0-9.,-]", "", table[4]), ",")[[1]]); 
# tables3$...6 = ifelse(tables3$...6 < 0.05,paste("<span style=\" font-weight: bold; \" >" ,sprintf("%.3f",tables3$...6), "</span>"),  paste("<span>" ,sprintf("%.3f",tables3$...6), "</span>"))
# colnames(tables3)[5] = "\u03C7\u00B2"; colnames(tables3)[6] = "p"
# tables3$p[1]= "<span style=\" font-weight: bold;    \" >\u003C 0.001</span>"
# tables3 %>%   
#   kbl(caption ="Latency (ms)",escape = F) %>%
#   kable_styling(latex_options = "HOLD_position", position = "center", full_width = F) %>%  row_spec(0,bold=T,align='c')
# tmp = tables2[(length(table$Effect)+1):(length(table$Effect)+5),1:2]
# names(tmp) <- NULL
# tmp1 <- data.frame(t(tmp[-1]))
# colnames(tmp1) <- tmp[[1]]
# tmp1 %>% kbl(digits = 2) %>%
#   kable_styling(latex_options = "HOLD_position", position = "center", full_width = F) 


# -------------------------------------- Regression table summary --------------------------------------

tab_model(bmod_full, show.p = F, show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 3, dv.labels = "Latency (ms)", file = "tmp/temp2.html", transform = NULL,  rm.terms = "beta")

report = capture.output(sexit(bmod_full, ci=.9))
print(paste("Bayesian general linear mixed model (exgaussian family with a identity link) (estimated using MCMC sampling with " ,chains ,"chains of", niter, "iterations and a warmup of", warm, ") to predict Latency with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: RT ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff). The model included condition, session, id and trialxcondition as random effects (formula: list(~condition * session | id, ~1 | trialxcondition)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 130.20) distributions"))
## [1] "Bayesian general linear mixed model (exgaussian family with a identity link) (estimated using MCMC sampling with  4 chains of 5000 iterations and a warmup of 1000 ) to predict Latency with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: RT ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff). The model included condition, session, id and trialxcondition as random effects (formula: list(~condition * session | id, ~1 | trialxcondition)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 130.20) distributions"
full_tab
## Summary of Posterior Distribution
## 
## Parameter                      | Median |  MAD |         90% CI |     pd |  Rhat |    BF
## ----------------------------------------------------------------------------------------
## condition                      |  -5.20 | 1.88 | [-8.19, -1.99] | 99.49% | 1.000 | 21.65
## intervention                   |   0.66 | 2.93 | [-4.12,  5.47] | 58.82% | 1.000 | 0.974
## session                        |   0.36 | 2.60 | [-3.91,  4.46] | 55.47% | 1.000 | 0.904
## age                            |   1.88 | 2.91 | [-2.77,  6.83] | 74.49% | 1.000 |  1.20
## gender                         |   0.38 | 2.94 | [-4.45,  5.05] | 54.83% | 1.000 | 0.987
## BMI_V1                         |   0.74 | 2.91 | [-3.98,  5.55] | 60.19% | 1.000 | 0.957
## hungry                         |  -2.10 | 2.70 | [-6.49,  2.41] | 78.28% | 1.000 |  1.14
## GLP_diff                       |  -0.90 | 2.91 | [-5.67,  3.96] | 61.81% | 1.000 |  1.04
## reelin_diff                    |  -0.19 | 2.98 | [-4.96,  4.83] | 52.59% | 1.000 | 0.985
## condition:intervention         |  -0.11 | 1.83 | [-3.11,  2.84] | 52.54% | 1.000 | 0.610
## condition:session              |  -2.02 | 1.72 | [-5.01,  0.82] | 87.66% | 1.000 |  1.20
## intervention:session           |  -1.30 | 2.50 | [-5.80,  2.51] | 70.23% | 1.000 | 0.956
## condition:intervention:session |   0.04 | 1.69 | [-2.83,  2.81] | 50.78% | 1.000 | 0.559
report[4]
## [1] "- b_condition (Median = -5.20, 90% CI [-8.19, -1.99]) has a 99.49% probability of being negative (< 0), 99.46% of being significant (< -0.05), and 99.24% of being large (< -0.30)"
tables <- list.clean(readHTMLTable("tmp/temp2.html"), fun = is.null, recursive = FALSE)
tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
tables2 <- as.matrix(tables2) %>% as_tibble()

# check for bad ICC?
tables2[is.na(tables2)] <- "";  names(tables2) <- NULL; tmp = t(tables2[(length(full_tab$Parameter)+1):(length(full_tab$Parameter)+5),1:2]);
tmp[,5][1] = "R2"; tmp[,5][2] = gsub(".*/","",tmp[,5][2]); colnames(tmp) =  tmp[1,]
pander:: pander(tmp[2,])
## 
## --------------------------------------------------------
##  ICC    N id   N trialxcondition   Observations    R2   
## ------ ------ ------------------- -------------- -------
##  1.00    50           20               3223       0.204 
## --------------------------------------------------------


Plot Latency

dfdraws = bmod_full %>%
    spread_draws(`b_condition` )


HDI_RT = plotHDI( dfdraws$`b_condition` , credMass = .90, binSize = 100, Title = "") + theme_bw() + html_theme  


dfdraws2 =  bmod_full %>%
     emmeans(~ condition) %>%
     gather_emmeans_draws() 


pp = dfdraws2 %>%
    ggplot(aes(x = as.factor(condition) , y = .value,  fill = as.factor(condition))) +
    #geom_abline(slope= 0, intercept=0, linetype = "dashed", color = "black") +
    stat_slab(.width = c(0.50, 0.9), position="dodge", alpha=0.5) +
    stat_pointinterval(.width = c(0.50, 0.9),position="dodge") +
    labs(x = "", y = "Latency (ms)", title = "") + 
    scale_fill_manual(values=c("-1" = pal[1],"1"=pal[2]), guide="none") +
    scale_color_manual( values=c("-1" = pal[1],"1"=pal[2]), guide="none") +
    scale_x_discrete(labels=c("CS-", "CS+")) +
  #scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(400,550, by = 50)), limits = c(380,575)) +
    theme_bw()

plt_bayes_html = pp + html_theme 
plt_bayes = pp + averaged_theme 


figureRT <- ggarrange(plt_bayes_html, HDI_RT,
                    labels = c("A", "B"),
                    ncol = 2)
#figureRT



# RT


dfR <- summarySEwithin(PAV.means,
                       measurevar = "RT",
                       withinvars = c("condition","session"), 
                       idvar = "id")

dfR$cond <- ifelse(dfR$condition == "1", -0.25, 0.25)
PAV.means$cond <- ifelse(PAV.means$condition == "1", -0.25, 0.25)
PAV.means <- PAV.means %>% mutate(condjit = jitter(as.numeric(cond), 0.3),
                                  grouping = interaction(id, cond))

pp <- ggplot(PAV.means, aes(x = cond, y = RT, 
                            fill = as.factor(condition), color = as.factor(condition))) +
  geom_line(aes(x = condjit, group = id, y = RT), alpha = .3, size = 0.5, color = 'gray') +
  geom_flat_violin(scale = "count", trim = FALSE, alpha = .2, aes(fill = as.factor(condition), color = NA))+
  geom_point(aes(x = condjit, shape = as.factor(intervention)), alpha = .3,) +
  geom_crossbar(data = dfR, aes(y = RT, ymin=RT-se, ymax=RT+se), width = 0.2 , alpha = 0.1)+
  ylab('Latency (ms)')+
  xlab('Conditioned stimulus')+
  scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(200,1000, by = 200)), limits = c(150,1050)) +
  scale_x_continuous(labels=c("CS+", "CS-"),breaks = c(-.25,.25), limits = c(-.5,.5)) +
  scale_fill_manual(values=c("1"= pal[2], "-1"=  pal[1]), guide = 'none') +
  scale_color_manual(values=c("1"= pal[2], "-1"=  pal[1]), guide = 'none') +
  scale_shape_manual(name="Intervention", labels=c("Placebo", "Liraglutide"), values = c(1, 2)) +
  theme_bw() + facet_wrap(~session, labeller=labeller(session = labels))

plt = pp + averaged_theme  + theme(legend.position=c(.96,.94), axis.text.x = element_text(size = 14))
pp + html_theme  + theme(legend.position=c(.96,.94), axis.text.x = element_text(size = 14))
A) Posterior distribution by Pavlovian cue. B) Highest density interval (90% HDI) of the posterior distribution <br> difference for the latency to respond between CS+ and CS-

  1. Posterior distribution by Pavlovian cue. B) Highest density interval (90% HDI) of the posterior distribution
    difference for the latency to respond between CS+ and CS- 

cairo_pdf('figures/Figure_PavlovianRT_Lira.pdf')
print(plt)
dev.off

cairo_pdf('figures/Figure_PavlovianRT_Lira_Bayes.pdf')
print(figureRT)
dev.off()

Perceived liking (Pavlovian Cue)

Ratings = how pleasant is the clue (0-100, no repetitions) & condition = CS+ or CS-

# Model
mf = formula(liking ~ condition*intervention*session + age + gender + BMI_V1  + hungry + GLP_diff + reelin_diff + (condition*session|id))

# -------------------------------------- Bayesian Regression Models using Stan -------------------------------------


bmod_full = brm(mf, data=PAV.means, family = gaussian, prior = c(prior(normal(0, 3), class = "b", coef = ""), prior(normal(0, 100), class = "Intercept", coef = "")), sample_prior = T, save_pars = save_pars(all = TRUE), chains = chains,  iter = niter, warmup = warm, seed = 123, inits = 1, backend = 'cmdstanr', threads = threading(4), control = list(adapt_delta = 0.99))  # a lot to unwind here..  1) Generic weakly informative prior around 0 for fixed effects and very weak prior for the intercept 2) we need to sample priors and save parameters for computing BF  3)larger step size  



#lmer to compare
fmod_full <- lmer(update(mf, ~ .-(condition*session|id)+(condition+session|id)) , data=PAV.means) #cannot maximize random structure here
## plot population-level effects posterior distributions and chain sampling

param = mcmc_plot(object = bmod_full, pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="areas")

trace = mcmc_plot(object = bmod_full, pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="trace")


#check assumptions
var_group = pp_check(bmod_full, type = "stat_grouped", group = "intervention", binwidth = 1, nsamples = NULL) #equality of variance between groups

rep_fit = pp_check(bmod_full, nsamples = 10) # check response fit 

error = pp_check(bmod_full, type ="error_scatter_avg", nsamples = NULL) # check good alignment between model and data, and no obvious pattern to the types of errors we are getting.


#Normality of errors
residuals <-residuals(bmod_full)[, 1]; res <- qqnorm(residuals, pch = 1, plot.it=F)

lmx = plot_model(fmod_full, type = "diag"); #diagnostic plots for lmer

#plot all
diagLik <- ggarrange(param, var_group, rep_fit, error, ncol = 2, nrow = 2)

annotate_figure(diagLik, top = text_grob("Diagnostic Plots", face = "bold", size = 10))

full_tab = describe_posterior(bmod_full,
                   estimate = "median", dispersion = T,
                   ci = .9, ci_method = "hdi",
                   bf_prior = bmod_full, diagnostic = "Rhat",  
                   test = c("p_direction", "bf"))
full_tab = filter(full_tab, Parameter %notin% c("b_Intercept", "prior_beta", "prior_sigma"))



# -------------------------------------- FREQUENTIST STATS: -----------------------------------------------

# measurevar = "liking"
# formula = "condition*intervention*session"
# COVA = "+ age + gender + BMI_V1 + "
# random = "Error(id/session*condition)"
#  
# mdl.lik = aov_car(as.formula(paste(measurevar, paste (formula, COVA, random), sep="~")), data= PAV.means, factorize = F, anova_table = list(correction = "GG", es = "pes"))
# res = nice(mdl.lik, MSE=F);  res = res[c(10,1:6,11,15,16),] #clean up unesuful an reorder
# ref_grid(mdl.lik)  #triple check everything is centered at 0
# 
# #calculate Partial eta-squared and its 90 % CI for each effect
# PES.lik = pes_ci(as.formula(paste(measurevar, paste (formula, COVA, random), sep="~")), data = PAV.means, conf.level = .90, factorize = FALSE, anova.type = "II", epsilon="none") ;  PES.lik = PES.lik[c(10,1:6,11,15,16),]
# mdl.lik.emms = emmeans(mdl.lik, pairwise ~ condition)
# 
# res$p.value = as.numeric(res$p.value)
# res$p.value = ifelse(res$p.value < 0.05,paste("<span style=\" font-weight: bold; \" >" ,sprintf("%.3f",res$p.value), "</span>"),  paste("<span>" ,sprintf("%.3f",res$p.value), "</span>"))
# 
# res$F = unlist(str_split(gsub("[^0-9.,-]", "", res$F), ","));res$pes = unlist(str_split(gsub("[^0-9.,-]", "", res$pes), ","));
# res$`90% CI` = paste(sprintf("%.3f",PES.lik[,2]), "-", sprintf("%.3f",PES.lik[,3]))
# 
# res$p.value[1]= "<span style=\" font-weight: bold;    \" >\u003C 0.001</span>"
# res$pes[c(5,7,8)]= c("\u003C 0.001")
# colnames(res)[3:5] = c( paste("F(", res$df[1], ")", sep=""),"&eta;<sub>p</sub><sup>2</sup>", "p")
# res[c(1,4,6,3,5)]  %>% kbl(digits = 2, escape = F,row.names = F)  %>%
#   kable_styling(latex_options = "striped", position = "center", full_width = F) 

# -------------------------------------- Regression table summary --------------------------------------

tab_model(fmod_full, show.p = F, show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 3, dv.labels = "Perceived liking", file = "tmp/temp3.html", rm.terms = "beta")

report = capture.output(sexit(bmod_full, ci=.9))
print(paste("Bayesian linear mixed model (estimated using MCMC sampling with" ,chains ," chains of", niter, " iterations and a warmup of", warm, ")  to predict liking with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: liking ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff). The model included condition, session and id as random effects (formula: ~condition * session | id). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 15.10) distributions"))
## [1] "Bayesian linear mixed model (estimated using MCMC sampling with 4  chains of 5000  iterations and a warmup of 1000 )  to predict liking with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: liking ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff). The model included condition, session and id as random effects (formula: ~condition * session | id). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 15.10) distributions"
full_tab
## Summary of Posterior Distribution
## 
## Parameter                      | Median |  MAD |         90% CI |     pd |  Rhat |     BF
## -----------------------------------------------------------------------------------------
## condition                      |   6.86 | 1.66 | [ 3.98,  9.48] | 99.99% | 1.000 | 235.33
## intervention                   |  -1.88 | 1.29 | [-4.01,  0.24] | 92.46% | 1.000 |   1.26
## session                        |  -1.94 | 1.03 | [-3.59, -0.21] | 97.03% | 1.000 |   2.12
## age                            |   1.52 | 1.22 | [-0.47,  3.52] | 89.61% | 1.000 |  0.898
## gender                         |   0.04 | 1.24 | [-2.03,  2.12] | 51.34% | 1.000 |  0.424
## BMI_V1                         |   0.38 | 1.29 | [-1.67,  2.56] | 61.96% | 1.000 |  0.459
## hungry                         |   4.64 | 1.16 | [ 2.78,  6.59] | 99.98% | 1.000 | 438.58
## GLP_diff                       |   0.38 | 1.62 | [-2.20,  3.06] | 59.21% | 1.000 |  0.551
## reelin_diff                    |  -0.39 | 1.64 | [-3.00,  2.39] | 59.71% | 1.000 |  0.549
## condition:intervention         |   0.39 | 1.57 | [-2.27,  2.93] | 59.68% | 1.000 |  0.542
## condition:session              |  -0.37 | 0.98 | [-2.04,  1.16] | 65.11% | 1.000 |  0.339
## intervention:session           |   0.32 | 1.00 | [-1.37,  1.99] | 62.52% | 1.000 |  0.342
## condition:intervention:session |   1.09 | 0.96 | [-0.44,  2.73] | 87.25% | 1.000 |  0.600
report[c(4,10)]
## [1] "- b_condition (Median = 6.86, 90% CI [3.98, 9.48]) has a 99.99% probability of being positive (> 0), 99.90% of being significant (> 1.05), and 63.54% of being large (> 6.29)"
## [2] "- b_hungry (Median = 4.64, 90% CI [2.78, 6.59]) has a 99.98% probability of being positive (> 0), 99.86% of being significant (> 1.05), and 7.46% of being large (> 6.29)"
tables <- list.clean(readHTMLTable("tmp/temp3.html"), fun = is.null, recursive = FALSE)
tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
tables2 <- as.matrix(tables2) %>% as_tibble()


tables2[is.na(tables2)] <- "";  names(tables2) <- NULL; tmp = t(tables2[(length(full_tab$Parameter)+1):(length(full_tab$Parameter)+4),1:2]);
tmp[,4][1] = "R2"; tmp[,4][2] = gsub(".*/","",tmp[,4][2])
pander::pander(tmp)
## 
## ------ ------ -------------- -------
##  ICC    N id   Observations    R2   
## 
##  0.55    50        184        0.681 
## ------ ------ -------------- -------



Plot Perceived liking (Pavlovian Cue)

dfdraws = bmod_full %>%
    spread_draws(`b_condition`,`b_hungry` )


HDI_cond_PAV = plotHDI( dfdraws$`b_condition` , credMass = .90, binSize = 100, Title = "") + theme_bw()

HDI_hunger_PAV = plotHDI( dfdraws$`b_hungry` , credMass = .90, binSize = 100, Title = "") + theme_bw() # also mcmc_plot(object = bmod_full, pars =c("b_hungry"), type ="areas", prob = 0.9)



dfdraws2 =  bmod_full %>%
     emmeans(~ condition) %>%
     gather_emmeans_draws() 


pp = dfdraws2 %>%
    ggplot(aes(x = as.factor(condition) , y = .value,  fill = as.factor(condition))) +
    geom_abline(slope= 0, intercept=0, linetype = "dashed", color = "black") +
    #geom_point(position = "jitter") +
    stat_slab(.width = c(0.50, 0.9), position="dodge", alpha=0.5) +
    stat_pointinterval(.width = c(0.50, 0.9),position="dodge") +
    ylab('Perceived liking')+
    xlab('')+
    scale_fill_manual(values=c("-1" = pal[1],"1"=pal[2]), guide="none") +
    scale_color_manual( values=c("-1" = pal[1],"1"=pal[2]), guide="none") +
    scale_x_discrete(labels=c("CS-", "CS+")) +
     scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(0,100, by = 25)), limits = c(-0.05,100.5)) +
    theme_bw()# + facet_wrap(~group, labeller=labeller(group =labels))

plt_bayes_html = pp + html_theme 
plt_bayes = pp + averaged_theme 

figurePAV <- ggarrange(plt_bayes_html, HDI_cond_PAV,
                    labels = c("A", "B"),
                    ncol = 2)
#figurePAV


# Liking

dfL <- summarySEwithin(PAV.means,
                       measurevar = "liking",
                       withinvars = c("condition", "session"), 
                       idvar = "id")

dfL$cond <- ifelse(dfL$condition == "1", -0.25, 0.25)


pp <- ggplot(PAV.means, aes(x = cond, y = liking, 
                            fill = as.factor(condition), color = as.factor(condition))) +
  geom_hline(yintercept=50, linetype="dashed", size=0.4, alpha=0.8) +
  geom_line(aes(x = condjit, group = id, y = liking), alpha = .3, size = 0.5, color = 'gray' ) +
  geom_flat_violin(scale = "count", trim = FALSE, alpha = .2, aes(fill =  as.factor(condition), color = NA)) +
  geom_point(aes(x = condjit, shape =  as.factor(intervention)), alpha = .3) +
  geom_crossbar(data = dfL, aes(y = liking, ymin=liking-se, ymax=liking+se), width = 0.2 , alpha = 0.1)+
  ylab('Liking Ratings')+
  xlab('Conditioned stimulus')+
  scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(0,100, by = 25)), limits = c(-0.05,100.5)) +
  scale_x_continuous(labels=c("CS+", "CS-"),breaks = c(-.25,.25), limits = c(-.5,.5)) +
  scale_fill_manual(values=c("1"= pal[2], "-1"=  pal[1]), guide = 'none') +
  scale_color_manual(values=c("1"= pal[2], "-1"=  pal[1]), guide = 'none') +
  scale_shape_manual(name="Intervention", labels=c("Placebo", "Liraglutide"), values = c(1, 2)) +
  theme_bw()+ facet_wrap(~session, labeller=labeller(session = labels))


plt = pp + averaged_theme  + theme(legend.position=c(.96,.94), axis.text.x = element_text(size = 14))
pp + html_theme  + theme(legend.position=c(.96,.94), axis.text.x = element_text(size = 14))
A) Posterior distribution by Pavlovian cue. B) Highest density interval (90% HDI) of the posterior distribution <br> difference for the latency to respond between CS+ and CS-

  1. Posterior distribution by Pavlovian cue. B) Highest density interval (90% HDI) of the posterior distribution
    difference for the latency to respond between CS+ and CS- 



cairo_pdf('figures/Figure_PavlovianLiking_Lira.pdf')
print(plt)
dev.off()

cairo_pdf('figures/Figure_PavlovianLiking_Lira_Bayes.pdf')
print(figurePAV)
dev.off()

Instrumental Conditioning Task

grips = number of times participant exceeded the force threshold to acquire the reward (Milkshake)

# Model
mf1 = formula(grips ~ trialZ*intervention*session + age + gender + BMI_V1  + hungry + GLP_diff + reelin_diff + (session|id) + (1|trialZ)) # LINEAR FIT  
mf2 = formula(grips ~  lspline(trialZ,-1.08327290)*intervention*session + age + gender + BMI_V1  + hungry + GLP_diff + reelin_diff + (session|id) + (1|trialZ)) # PIECEWISE REGRESSION WITH SPLINE AT 5 #(-1.08327290 is 5 after scaling)


# -------------------------------------- Bayesian Regression Models using Stan -------------------------------------
linmod = brm(mf1, data=INST, family = gaussian, prior = c(prior(normal(0, 3), class = "b", coef = ""), prior(normal(0, 100), class = "Intercept", coef = "")), sample_prior = T, save_pars = save_pars(all = TRUE), chains = chains,  iter = niter, warmup = warm, seed = 123, inits = 1, backend = 'cmdstanr', threads = threading(4))  # a lot to unwind here..  1) Generic weakly informative prior around 0 for fixed effects and very weak prior for the intercept 2) we need to sample priors and save parameters for computing BF  
splinemod = update(linmod, formula. = mf2)
linmod$loo = loo(linmod); splinemod$loo = loo(splinemod)
comp = loo::loo_compare(linmod$loo, splinemod$loo) 

bmod_full = splinemod #winner model

#lmer to compare
fmod_full = lmer(mf2 , data = INST, REML=F, control = control )

Model Comparison between linear and piecewise with spline regression
Best fit => Piecewise Regression with spline: smaller ELPD (expected log pointwise predictive density estimated
via leave-one-out cross-validation) is better

comp = tibble::rownames_to_column(as_tibble(comp, rownames = NA)); colnames(comp)[1] = "model"
pander(comp[c(1:3,8:9)]) # so the splinemod was preferred as the other model had smaller ELPD. Thus, spline improved out-of-sample predictions. can also see: compare_ic(linmod, splinemod) but its deprecated
## 
## ----------------------------------------------------
##    model     elpd_diff   se_diff   looic   se_looic 
## ----------- ----------- --------- ------- ----------
##  splinemod       0          0      11920    93.67   
## 
##   linmod      -1.373      2.866    11923    93.43   
## ----------------------------------------------------
## plot population-level effects posterior distributions and chain sampling

param = mcmc_plot(object = bmod_full, pars =c("b_.*al"), type ="areas")
param2 = mcmc_plot(object = bmod_full, pars =c("b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="areas")
param = param + scale_y_discrete(labels=c("trial<5", "trial>5", "trial<5:intervention", "trial>5:intervention","trial<5:session", "trial>5:session", "trial<5:intervention:session", "trial>5:intervention:session"))

trace = mcmc_plot(object = bmod_full, pars =c("b_.*al"), type ="trace")
trace2 = mcmc_plot(object = bmod_full, pars =c("b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="trace")

#check assumptions
var_group = pp_check(bmod_full, type = "stat_grouped", group = "intervention", binwidth = 0.01, nsamples = NULL) #equality of variance between groups

rep_fit = pp_check(bmod_full, nsamples = 10) # check response fit 

error = pp_check(bmod_full, type ="error_scatter_avg", nsamples = NULL) # check good alignment between model and data, and no obvious pattern to the types of errors we are getting.


#Normality of errors
residuals <-residuals(bmod_full)[, 1]; res <- qqnorm(residuals, pch = 1, plot.it=F)

lmx = plot_model(fmod_full, type = "diag"); #diagnostic plots for lmer

#plot all
diagINST <- ggarrange(param, var_group, rep_fit, error, ncol = 2, nrow = 2)

annotate_figure(diagINST, top = text_grob("Diagnostic Plots", face = "bold", size = 10))

full_tab = describe_posterior(bmod_full,
                   estimate = "median", dispersion = T,
                   ci = .9, ci_method = "hdi",
                   bf_prior = bmod_full, diagnostic = "Rhat",  
                   test = c("p_direction", "bf"))
full_tab = filter(full_tab, Parameter %notin% c("b_Intercept", "prior_beta", "prior_sigma"))

x = attributes(full_tab); x$clean_parameters[c(2:3,12:15,17:18),5] = c("trial<5",       "trial>5", "trial<5:intervention", "trial>5:intervention","trial<5:session", "trial>5:session", "trial<5:intervention:session", "trial>5:intervention:session"); attributes(full_tab) = x


# -------------------------------------- FREQUENTIST STATS: LRT + Bootstrap  -----------------------------------------------

# fmod_spline =  update(fmod_full, ~ .-lspline(trial,-1.08327290)) #evaluate trial spline
# 
# # p-value from bootstrap distribution
# LRT_spline = PBmodcomp(fmod_full, fmod_spline, nsim = 500, seed = 123, cl=cores)


# -------------------------------------- Regression table summary --------------------------------------

tab_model(bmod_full, show.p = F, show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 3, dv.labels = "Number of Grips", file = "tmp/temp4.html", transform = NULL,  rm.terms = "beta")

report = capture.output(sexit(bmod_full, ci=.9, large = 0.69))
print(paste("Bayesian linear mixed model (estimated using MCMC sampling with" ,chains ," chains of", niter, " iterations and a warmup of", warm, ") to predict grips with trial, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: grips ~ lspline(trial, 5) + intervention + session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff + lspline(trial, 5):intervention + lspline(trial, 5):session + intervention:session + lspline(trial, 5):intervention:session). The model included session, id and trial as random effects (formula: list(~session | id, ~1 | trial)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 8.90) distributions"))
## [1] "Bayesian linear mixed model (estimated using MCMC sampling with 4  chains of 5000  iterations and a warmup of 1000 ) to predict grips with trial, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: grips ~ lspline(trial, 5) + intervention + session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff + lspline(trial, 5):intervention + lspline(trial, 5):session + intervention:session + lspline(trial, 5):intervention:session). The model included session, id and trial as random effects (formula: list(~session | id, ~1 | trial)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 8.90) distributions"
full_tab
## Summary of Posterior Distribution
## 
## Parameter                    | Median |  MAD |         90% CI |     pd |  Rhat |     BF
## ---------------------------------------------------------------------------------------
## trial<5                      |   0.39 | 0.61 | [-0.65,  1.38] | 73.96% | 1.000 |  0.256
## trial>5                      |  -0.82 | 0.10 | [-0.99, -0.65] |   100% | 1.000 | > 1000
## intervention                 |  -0.58 | 1.12 | [-2.47,  1.22] | 70.38% | 1.000 |  0.430
## session                      |  -0.95 | 0.78 | [-2.22,  0.33] | 88.96% | 1.000 |  0.564
## age                          |  -1.07 | 0.85 | [-2.48,  0.33] | 89.60% | 1.000 |  0.635
## gender                       |  -0.27 | 0.88 | [-1.73,  1.19] | 61.84% | 1.000 |  0.300
## BMI_V1                       |  -1.37 | 0.95 | [-2.96,  0.22] | 91.88% | 1.000 |  0.894
## hungry                       |  -0.17 | 0.68 | [-1.29,  0.97] | 59.91% | 1.000 |  0.233
## GLP_diff                     |   1.11 | 1.19 | [-0.99,  3.00] | 82.27% | 1.000 |  0.618
## reelin_diff                  |   0.21 | 1.31 | [-1.88,  2.41] | 56.27% | 1.000 |  0.442
## trial<5:intervention         |   0.69 | 0.54 | [-0.20,  1.56] | 89.87% | 1.000 |  0.397
## trial>5:intervention         |   0.17 | 0.09 | [ 0.03,  0.32] | 96.84% | 1.000 |  0.181
## trial<5:session              |  -1.26 | 0.54 | [-2.12, -0.37] | 98.94% | 1.000 |   2.93
## trial>5:session              |   0.13 | 0.09 | [-0.02,  0.28] | 91.89% | 1.000 |  0.078
## intervention:session         |  -1.15 | 0.78 | [-2.41,  0.15] | 92.96% | 1.000 |  0.764
## trial<5:intervention:session |  -0.75 | 0.53 | [-1.65,  0.11] | 92.34% | 1.000 |  0.505
## trial>5:intervention:session |  -0.17 | 0.09 | [-0.32, -0.03] | 97.43% | 1.000 |  0.204
x = report[c(4:5,16)]; x = gsub('b_lsplinetrialZM1.08327291', 'b_trial<5', x); x = gsub('b_lsplinetrialZM1.08327292', 'b_trial>5', x); x
## [1] "- b_trial<5 (Median = 0.39, 90% CI [-0.65, 1.38]) has a 73.96% probability of being positive (> 0), 51.01% of being significant (> 0.37), and 31.08% of being large (> 0.69)"            
## [2] "- b_trial>5 (Median = -0.82, 90% CI [-0.99, -0.65]) has a 100.00% probability of being negative (< 0), 100.00% of being significant (< -0.37), and 89.06% of being large (< -0.69)"      
## [3] "- b_trial<5.session (Median = -1.26, 90% CI [-2.12, -0.37]) has a 98.94% probability of being negative (< 0), 94.95% of being significant (< -0.37), and 84.99% of being large (< -0.69)"
tables <- list.clean(readHTMLTable("tmp/temp4.html"), fun = is.null, recursive = FALSE)
tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
tables2 <- as.matrix(tables2) %>% as_tibble()


tables2[is.na(tables2)] <- "";  names(tables2) <- NULL; tmp = t( tables2[(length(full_tab$Parameter)+1):(length(full_tab$Parameter)+5),1:2]);
tmp[,5][1] = "R2"; tmp[,5][2] = gsub(".*/","",tmp[,5][2])
pander::pander(tmp)
## 
## ------ ------ ---------- -------------- -------
##  ICC    N id   N trialZ   Observations    R2   
## 
##  0.79    50       24          2232       0.789 
## ------ ------ ---------- -------------- -------



Plot

int_conditions <- list(
  session = setNames(c(-1, 1), c("pre", "post"))
)
plt = conditional_effects(bmod_full, effects = "trialZ:session",
                    int_conditions = int_conditions) #easier to extract spline with brms

dfTRIAL = as_tibble(plt$trial); dfTRIAL$grips = dfTRIAL$estimate__ ; dfTRIAL$trial = rescale(dfTRIAL$trialZ, to = c(1,24))


#by groups

dfTRIALg <- summarySEwithin(INST.means,
                            measurevar = "grips",
                            withinvars = c("trial", "session"),
                            betweenvars = "intervention",
                            idvar = "id")

dfTRIALg$trial       <- as.numeric(as.character(dfTRIALg$trial))

pp <-  ggplot(dfTRIAL, aes(x = trial, y = grips)) +
  geom_point(data = dfTRIALg, aes(shape = intervention), alpha = 0.3, color = 'black') +
  geom_line(data = dfTRIAL, size =1, color = pal[4]) +
  geom_ribbon(aes(ymin=grips-se__, ymax=grips+se__), fill = pal[4], alpha = 0.3, )+
  ylab('Number of Grips')+
  xlab('Trial') +
  scale_y_continuous(expand = c(0, 0),  limits = c(9.5,16.05),  breaks=c(seq.int(9,16, by = 1))) +
  #scale_x_continuous(expand = c(0, 0),  limits = c(-2,2),  breaks=c(seq.int(-2,2, by = 0.5))) +
  scale_shape_manual(name="Group", labels=c("Placebo", "Liraglutide"), values = c(1, 2, 18)) +
  theme_bw() +
  facet_wrap(~session, labeller=labeller(session = labels))



pp + html_theme + theme(legend.position=c(.9,.88),axis.text.x = element_text(size = 14))
Number of grips over trials by session.

Number of grips over trials by session. 

plt = pp + averaged_theme + theme(legend.position=c(.9,.88),axis.text.x = element_text(size = 14))



cairo_pdf('figures/Figure_INST_trial_LIRA.pdf')
print(plt)
dev.off()

Pavlvovian-Instrumental Transfer (PIT) Task

Mobilized effort = Area under the curve of the force exerted exceeding the delivery threshold during Pavlovian cue presentation

mf = formula(AUC ~ condition*intervention*session + age + gender + BMI_V1  + hungry + GLP_diff + reelin_diff + (condition*session|id) + (1|trialxcondition))

# -------------------------------------- Bayesian Regression Models using Stan -------------------------------------


bmod_full = brm(bf(mf, hu ~ 1), family = hurdle_gaussian, stanvars = stanvars, data=PIT.clean, prior =  c(prior(normal(0, 3), class = "b", coef = ""), prior(normal(0, 100), class = "Intercept", coef = ""), prior(logistic(0, 0.5), class = "Intercept", dpar = "hu")), sample_prior = T, save_pars = save_pars(all = TRUE), chains = chains,  iter = niter, warmup = warm, seed = 123, inits = 1, backend = 'cmdstan', threads = threading(4), control = list(adapt_delta = 0.99)) # a lot to unwind here.. 1) custom gaussian hurdle cause zero-inflated continous data more details here: https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html 2) Generic weakly informative prior around 0 for fixed effects and very weak prior for the intercept 3) we need to sample priors and save parameters for computing BF 4)increased step size

#lmer to compare
fmod_full = lmer(mf , data = PIT.clean, REML=F, control = control)
## plot population-level effects posterior distributions and chain sampling
param = mcmc_plot(object = bmod_full,  pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="areas") 

trace = mcmc_plot(object = bmod_full,  pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1"), type ="trace")


#check assumptions
var_group = pp_check(bmod_full, type = "stat_grouped", stat = "median", group = "intervention", binwidth = 1, nsamples = NULL) #equality of variance between groups

rep_fit = pp_check(bmod_full, nsamples = 10) # check response fit 

error = pp_check(bmod_full, type ="error_scatter_avg", nsamples = NULL) # check good alignment between model and data, and no obvious pattern to the types of errors we are getting.


#Normality of errors
residuals <-residuals(bmod_full)[, 1]; res <- qqnorm(residuals, pch = 1, plot.it=F)

lmx = plot_model(fmod_full, type = "diag"); #diagnostic plots for lmer

#plot all
diagPIT <- ggarrange(param, var_group, rep_fit, error, ncol = 2, nrow = 2)

annotate_figure(diagPIT, top = text_grob("Diagnostic Plots", face = "bold", size = 10))

# check response fit -> capture better the bimodal distrib
# residual ... could be better but at least it's MUCH better than simple gaussian see lmx[1] ... this is just catastrophic..
full_tab = describe_posterior(bmod_full,
                   estimate = "median", dispersion = T,
                   ci = .9, ci_method = "hdi",
                   bf_prior = bmod_full, diagnostic = "Rhat",  
                   test = c("p_direction", "bf"))
full_tab = filter(full_tab, Parameter %notin% c("b_Intercept", "b_hu_Intercept", "prior_sigma"))

# -------------------------------------- FREQUENTIST STATS: LRT + Bootstrap  -----------------------------------------------

 fmod_cond = update(fmod_full,  ~ .-condition) #evaluate condition

# # p-value from bootstrap distribution
 LRT_cond = PBmodcomp(fmod_full, fmod_cond, nsim = 500, seed = 123, cl=cores)

# -------------------------------------- Regression table summary --------------------------------------

tab_model(bmod_full, show.p = F,show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 3, dv.labels = "Mobilized effort (a.u.)", file = "tmp/temp5.html", transform = NULL,  rm.terms = "hu_Intercept")

report = capture.output(sexit(bmod_full, ci=.9))
print(paste("Bayesian general linear mixed model (hurdle gaussian family with a identity link) (estimated using MCMC sampling with" ,chains ," chains of", niter, " iterations and a warmup of", warm, ") to predict Mobilized Effort (AUC) with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: AUC ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff). The model included condition, session, id and trialxcondition as random effects (formula: list(~condition * session | id, ~1 | trialxcondition)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 130.20) distributions"))
## [1] "Bayesian general linear mixed model (hurdle gaussian family with a identity link) (estimated using MCMC sampling with 4  chains of 5000  iterations and a warmup of 1000 ) to predict Mobilized Effort (AUC) with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: AUC ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff). The model included condition, session, id and trialxcondition as random effects (formula: list(~condition * session | id, ~1 | trialxcondition)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00), and student_t (location = 0.00, scale = 130.20) distributions"
full_tab
## Summary of Posterior Distribution
## 
## Parameter                      | Median |  MAD |        90% CI |     pd |  Rhat |    BF
## ---------------------------------------------------------------------------------------
## condition                      |   3.12 | 2.10 | [-0.36, 6.54] | 93.10% | 1.001 |  2.05
## intervention                   |   0.62 | 2.93 | [-4.26, 5.39] | 58.22% | 1.000 | 0.992
## session                        |   0.06 | 2.75 | [-4.52, 4.63] | 50.79% | 1.000 | 0.925
## age                            |   1.68 | 2.88 | [-3.17, 6.37] | 71.89% | 1.000 |  1.13
## gender                         |   0.30 | 2.88 | [-4.44, 5.09] | 53.81% | 1.000 | 0.976
## BMI_V1                         |   0.06 | 2.87 | [-4.68, 4.89] | 50.96% | 1.000 | 0.937
## hungry                         |  -0.14 | 2.83 | [-4.83, 4.54] | 51.96% | 1.000 | 0.953
## GLP_diff                       |   0.89 | 2.97 | [-3.87, 5.90] | 61.94% | 1.000 |  1.04
## reelin_diff                    |   0.44 | 2.96 | [-4.47, 5.28] | 55.67% | 1.000 | 0.998
## condition:intervention         |   0.72 | 2.10 | [-2.69, 4.19] | 63.44% | 1.001 | 0.758
## condition:session              |  -0.11 | 1.73 | [-2.96, 2.72] | 52.58% | 1.000 | 0.570
## intervention:session           |  -0.81 | 2.72 | [-5.20, 3.81] | 61.92% | 1.000 | 0.925
## condition:intervention:session |   0.76 | 1.71 | [-2.09, 3.53] | 66.88% | 1.000 | 0.616
report[5]
## [1] "- b_condition (Median = 3.12, 90% CI [-0.36, 6.54]) has a 93.10% probability of being positive (> 0), 92.80% of being significant (> 0.05), and 91.04% of being large (> 0.30)"
tables <- list.clean(readHTMLTable("tmp/temp5.html"), fun = is.null, recursive = FALSE)
tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
tables2 <- as.matrix(tables2) %>% as_tibble(); 


tables2[is.na(tables2)] <- "";  names(tables2) <- NULL; tmp = t(tables2[(length(full_tab$Parameter)+1):(length(full_tab$Parameter)+5),1:2]);
tmp[,5][1] = "R2"; tmp[,5][2] = gsub(".*/","",tmp[,5][2])
pander::pander(tmp)
## 
## ------ ------ ------------------- -------------- -------
##  ICC    N id   N trialxcondition   Observations    R2   
## 
##  0.70    50           15               2760       0.709 
## ------ ------ ------------------- -------------- -------

Plot

# plot HDI
dfdraws = bmod_full %>%
    spread_draws(`b_condition` )

HDI_PIT = plotHDI( dfdraws$`b_condition` , credMass = .90, binSize = 100, Title = "") + theme_bw()


#plot estimated means from posterior distribution from the model draws


dfdraws2 =  bmod_full %>%
     emmeans(~ condition) %>%
     gather_emmeans_draws() 

# CSp = subset(dfdraws2, condition ==1); CSm = subset(dfdraws2, condition ==-1); diff= CSp; diff$.value = CSp$.value - CSm$.value


pp = dfdraws2 %>% #diff
    ggplot(aes(x = as.factor(condition) , y = .value,  fill = as.factor(condition))) +
    #geom_abline(slope= 0, intercept=0, linetype = "dashed", color = "black") +
    #geom_point(position = "jitter") +
    stat_slab(.width = c(0.50, 0.9), position="dodge", alpha=0.5) +
    stat_pointinterval(.width = c(0.50, 0.9),position="dodge") +
    ylab('Mobilized effort')+
    xlab('')+
    scale_fill_manual(values=c("1" = pal[6],"-1"=pal[1]), guide="none") +
    scale_color_manual( values=c("1" = pal[6],"-1"=pal[1]), guide="none") +
    scale_x_discrete(labels=c("CS-", "CS+")) +
    theme_bw()

plt_bayes_html = pp + html_theme 
plt_bayes = pp + averaged_theme 


figurePIT <- ggarrange(plt_bayes_html, HDI_PIT,
                    labels = c("A", "B"),
                    ncol = 2)
#figurePIT


#### Plot PIT

dfL <- summarySEwithin(PIT.means,
                       measurevar = "AUC",
                       withinvars = c("condition","session"), 
                       idvar = "id")

dfL$cond <- ifelse(dfL$condition == "1", -0.25, 0.25)
PIT.means$cond <- ifelse(PIT.means$condition == "1", -0.25, 0.25)
PIT.means <- PIT.means %>% mutate(condjit = jitter(as.numeric(cond), 0.3),grouping = interaction(id, cond))


# AUC
pp <- ggplot(PIT.means, aes(x = cond, y = AUC,
                            fill = as.factor(condition), color = as.factor(condition))) +
  geom_line(aes(x = condjit, group = id, y = AUC), alpha = .3, size = 0.5, color = 'gray' ) +
  geom_flat_violin(scale = "count", trim = FALSE, alpha = .2, aes(fill = as.factor(condition), color = NA)) +
  geom_point(aes(x = condjit, shape = as.factor(intervention)), alpha = .3) +
  geom_crossbar(data = dfL, aes(y = AUC, ymin=AUC-se, ymax=AUC+se), width = 0.2 , alpha = 0.1)+
  ylab('Mobilized effort (a.u.)')+
  xlab('Condition')+
  #scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(0,30, by = 5)), limits = c(-0.05,30.5)) +
  scale_x_continuous(labels=c("CS+", "CS-"),breaks = c(-.25,.25), limits = c(-.5,.5)) +
  scale_fill_manual(values=c("1"= pal[6], "-1"=  pal[1]), guide = 'none') +
  scale_color_manual(values=c("1"= pal[6], "-1"=  pal[1]), guide = 'none') +
  scale_shape_manual(name="Intervention", labels=c("Placebo", "Liraglutide"), values = c(1, 2)) +
  theme_bw()+ facet_wrap(~session, labeller=labeller(session = labels))



pp + html_theme +theme(legend.position=c(.94,.94))
A) Posterior distribution by group. B) Highest density interval (90% HDI) for the interaction condition by session

  1. Posterior distribution by group. B) Highest density interval (90% HDI) for the interaction condition by session 

plt = pp + averaged_theme +theme(legend.position=c(.94,.94))
# PLOT OVERTIME


pp <- ggplot(PIT.p, aes(x = as.numeric(trial), y = AUC,
                        color = condition,
                        fill  = condition))+
    geom_point(data = PIT.group, aes(shape = intervention, color = condition), alpha = 0.3) +
  geom_line(alpha = .5, size = 1, show.legend = F) +
  geom_ribbon(aes(ymax = AUC + se, ymin = AUC - se),  alpha=0.4) +
  geom_point() +
  ylab('Mobilized effort (a.u.)')+
  xlab('Trial')+
  scale_color_manual(labels = c('-1'= 'CS-', "1" = 'CS+'), name="",
                     values = c("1"= pal[6], '-1'= pal[1])) +
  scale_fill_manual(labels = c('-1'= 'CS-', "1" = 'CS+'), name="",
                    values = c("1"= pal[6], '-1'= pal[1])) +
  scale_y_continuous(expand = c(0, 0),  limits = c(50,200),  breaks=c(seq.int(50,200, by = 50))) +
  #scale_x_continuous(expand = c(0, 0),  limits = c(0,15),  breaks=c(seq.int(1,15, by = 2))) +
    scale_shape_manual(name="Group", labels=c("Placebo", "Liraglutide"), values = c(1, 2, 18)) +
  theme_bw() +
  facet_wrap(~session, labeller=labeller(session =labels))


pp + html_theme + theme(strip.background = element_rect(fill="white"), legend.key.size = unit(0.3, "cm"))
Mobilized effort for each condition and group over trials.

Mobilized effort for each condition and group over trials. 

plt = pp + averaged_theme + theme(strip.background = element_rect(fill="white"), legend.key.size = unit(0.8, "cm"), axis.text.x = element_text(size = 16))



# bmod_time = update(bmod_full,  ~ .+as.factor(trialxcondition)) #evaluate interaction
# 
# df =  bmod_time %>%
#      emmeans(~ condition:session:trialxcondition) 
# 
# 
# 
# bayes_plt = as_tibble(df) %>%
#     ggplot(aes(x = trialxcondition, y = emmean,  fill =as.factor(condition))) +
#     #geom_point(position = "jitter") +
#     geom_smooth(method=loess, se = T,inherit.aes = T) + 
#     #stat_slab(.width = c(0.50, 0.95),position="dodge", alpha=0.5) +
#     #stat_pointinterval(.width = c(0.50, 0.95),position="dodge") +
#     ylab('Mobilized effort (a.u.) \n Baseline corrected')+
#     xlab('')+
#     scale_fill_manual(values=c("1" = pal[2],"-1"=pal[1]), guide="none") +
#     scale_color_manual( values=c("1" = pal[2],"-1"=pal[1]), guide="none") +
#     #scale_x_discrete(labels=c("CS+", "CS-")) +
#     theme_bw() + facet_wrap(~session, labeller=labeller(session =labels))


cairo_pdf('figures/Figure_PIT_Lira.pdf')
print(plt)
dev.off()

cairo_pdf('figures/Figure_PIT_Lira_Bayes.pdf')
print(figurePIT)
dev.off()

Hedonic Reactivity Test

Perceived liking = how pleasant is the liquid solution rated (0-100, with repetitions)
&condition = Milshake or Tasteless & intensity = difference on how intense the liquid solution
were rated (mean(Milshake) - mean(Tasteless)) & familiarity = difference on how familiar the liquid solution
were rated (mean(Milshake) - mean(Tasteless))

mf = formula(perceived_liking ~ condition*intervention*session + age + gender + BMI_V1  + hungry + GLP_diff + reelin_diff + int + fam + (condition*session|id) + (1|trialxcondition))

# -------------------------------------- Bayesian Regression Models using Stan -------------------------------------


bmod_full = brm(bf(mf, hu ~ 1), family = hurdle_gaussian, stanvars = stanvars, data=HED, prior =  c(prior(normal(0, 3), class = "b", coef = ""), prior(normal(0, 100), class = "Intercept", coef = ""), prior(logistic(0, 0.5), class = "Intercept", dpar = "hu")), sample_prior = T, save_pars = save_pars(all = TRUE), chains = chains,  iter = niter, warmup = warm, seed = 123, inits = 1, backend = 'cmdstan', threads = threading(4), control = list(adapt_delta = 0.99))# a lot to unwind here.. 1)  custom gaussian hurdle cause zero-inflated continous data 2) Generic weakly informative prior around 0 for fixed effects and very weak prior for the intercept 3) we need to sample priors and save parameters for computing BF #this one is a big longer, so seat tight 4) increased step size
#lmer to compare
fmod_full = lmer(mf , data = HED, REML=F, control = control)
## plot population-level effects posterior distributions and chain sampling
param = mcmc_plot(object = bmod_full,  pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1", "b_int", "b_fam"), type ="areas")  # int and fam

trace = mcmc_plot(object = bmod_full,  pars =c("b_.*o", "b_.*y", "b_.*ge", "b_.*diff", "b_.*V1", "b_int", "b_fam"), type ="trace")


#check assumptions
var_group = pp_check(bmod_full, type = "stat_grouped", stat = "median", group = "intervention", binwidth = 0.1, nsamples = NULL) #equality of variance between groups

rep_fit = pp_check(bmod_full, nsamples = 10) # check response fit 

error = pp_check(bmod_full, type ="error_scatter_avg", nsamples = NULL) # check good alignment between model and data, and no obvious pattern to the types of errors we are getting.


#Normality of errors
residuals <-residuals(bmod_full)[, 1]; res <- qqnorm(residuals, pch = 1, plot.it=F)

lmx = plot_model(fmod_full, type = "diag"); #diagnostic plots for lmer

#plot all
diagPIT <- ggarrange(param, var_group, rep_fit, error, ncol = 2, nrow = 2)

annotate_figure(diagPIT, top = text_grob("Diagnostic Plots", face = "bold", size = 10))

# check response fit -> capture better the weird (bimodal at least..) distrib
# residual ... could be better but at least it's MUCH better than simple gaussian see lmx[1] ... this is just catastrophic..
full_tab = describe_posterior(bmod_full,
                   estimate = "median", dispersion = T,
                   ci = .9, ci_method = "hdi",
                   bf_prior = bmod_full, diagnostic = "Rhat",  
                   test = c("p_direction", "bf"))
full_tab = filter(full_tab, Parameter %notin% c("b_Intercept", "b_hu_Intercept", "prior_sigma"))


#contrasts 
# emmip(bmod_full, condition~intervention) to visualize
# 
# ems = emmeans(bmod_full, ~condition|intervention)
# con_inter = contrast(ems, interaction = "pairwise", by = NULL)
# 
# 
# inter_tab = describe_posterior(con_inter,
#                    estimate = "median", dispersion = TRUE,
#                    ci = .9, ci_method = "hdi",
#                    bf_prior = bmod_full, diagnostic = "Rhat",  
#                    test = c("p_direction", "bf"))



# -------------------------------------- FREQUENTIST STATS: LRT + Bootstrap  -----------------------------------------------

# fmod_cond = update(fmod_full,  ~ .-condition) #evaluate condition

# # p-value from bootstrap distribution
# LRT_cond = PBmodcomp(fmod_full, fmod_cond, nsim = 500, seed = 123, cl=cores) 

# -------------------------------------- Regression table summary --------------------------------------

tab_model(bmod_full, show.p = F,show.intercept = F, show.se = T, title ="", show.re.var = F, digits = 3, dv.labels = "Perceived liking", file = "tmp/temp6.html", transform = NULL,  rm.terms = "hu_Intercept")

report = capture.output(sexit(bmod_full, ci=.9))
print(paste("Bayesian general linear mixed model (hurdle gaussian family with a identity link) (estimated using MCMC sampling with", chains," chains of", niter, " iterations and a warmup of", warm, ") to predict Perceived Liking with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: perceived_liking ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff + int + fam). The model included condition, session, id and trialxcondition as random effects (formula: list(~condition | id, ~1 | trialxcondition)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00) and student_t (location = 0.00, scale = 31.7) distributions for beta and sd respectively"))
## [1] "Bayesian general linear mixed model (hurdle gaussian family with a identity link) (estimated using MCMC sampling with 4  chains of 5000  iterations and a warmup of 1000 ) to predict Perceived Liking with condition, intervention, session, age, gender, BMI_V1, hungry, GLP_diff and reelin_diff (formula: perceived_liking ~ condition * intervention * session + age + gender + BMI_V1 + hungry + GLP_diff + reelin_diff + int + fam). The model included condition, session, id and trialxcondition as random effects (formula: list(~condition | id, ~1 | trialxcondition)). Priors over parameters were set as normal (mean = 0.00, SD = 3.00) and student_t (location = 0.00, scale = 31.7) distributions for beta and sd respectively"
full_tab
## Summary of Posterior Distribution
## 
## Parameter                      | Median |  MAD |         90% CI |     pd |  Rhat |     BF
## -----------------------------------------------------------------------------------------
## condition                      |  12.91 | 1.75 | [10.11, 15.79] |   100% | 1.001 | > 1000
## intervention                   |   1.36 | 1.68 | [-1.52,  4.08] | 78.32% | 1.000 |  0.771
## session                        |   1.23 | 0.96 | [-0.33,  2.86] | 89.62% | 1.000 |  0.734
## age                            |  -0.51 | 1.61 | [-3.14,  2.16] | 62.71% | 1.000 |  0.568
## gender                         |  -2.24 | 1.62 | [-4.96,  0.44] | 91.09% | 1.000 |   1.42
## BMI_V1                         |   2.05 | 1.68 | [-0.83,  4.77] | 88.71% | 1.000 |   1.18
## hungry                         |   1.36 | 1.31 | [-0.84,  3.50] | 84.78% | 1.000 |  0.743
## GLP_diff                       |  -0.58 | 1.91 | [-3.71,  2.56] | 62.38% | 1.000 |  0.623
## reelin_diff                    |  -2.25 | 2.10 | [-5.73,  1.23] | 85.66% | 1.000 |   1.25
## int                            |   0.86 | 2.57 | [-3.47,  5.03] | 63.09% | 1.001 |  0.909
## fam                            |   0.81 | 1.57 | [-1.76,  3.40] | 69.80% | 1.000 |  0.604
## condition:intervention         |   2.44 | 1.53 | [-0.15,  4.88] | 94.21% | 1.001 |   1.94
## condition:session              |  -0.64 | 0.90 | [-2.06,  0.97] | 76.19% | 1.001 |  0.388
## intervention:session           |   0.15 | 0.90 | [-1.30,  1.66] | 56.86% | 1.000 |  0.305
## condition:intervention:session |   0.57 | 0.85 | [-0.87,  1.98] | 75.19% | 1.000 |  0.363
report[5]
## [1] "- b_condition (Median = 12.91, 90% CI [10.11, 15.79]) has a 100.00% probability of being positive (> 0), 100.00% of being significant (> 0.05), and 100.00% of being large (> 0.30)"
tables <- list.clean(readHTMLTable("tmp/temp6.html"), fun = is.null, recursive = FALSE)
tables2 = tables[[1]] %>% janitor::row_to_names(row_number = 1)
tables2 <- as.matrix(tables2) %>% as_tibble(); 

tables2[is.na(tables2)] <- "";  names(tables2) <- NULL; tmp = t(tables2[(length(full_tab$Parameter)+1):(length(full_tab$Parameter)+5),1:2]);
tmp[,5][1] = "R2"; tmp[,5][2] = gsub(".*/","",tmp[,5][2])
pander::pander(tmp)
## 
## ------ ------ ------------------- -------------- -------
##  ICC    N id   N trialxcondition   Observations    R2   
## 
##  0.47    50           20               3680       0.798 
## ------ ------ ------------------- -------------- -------

Plot

# plot HDI
dfdraws = bmod_full %>%
    spread_draws(`b_condition` )

HDI_HED = plotHDI( dfdraws$`b_condition` , credMass = .90, binSize = 100, Title = "") + theme_bw()


#plot estimated means from posterior distribution from the model draws


dfdraws2 =  bmod_full %>%
     emmeans(~ condition) %>%
     gather_emmeans_draws() 


pp = dfdraws2 %>% #diff
    ggplot(aes(x = as.factor(condition) , y = .value,  fill = as.factor(condition))) +
    #geom_abline(slope= 0, intercept=0, linetype = "dashed", color = "black") +
    #geom_point(position = "jitter") +
    stat_slab(.width = c(0.50, 0.9), position="dodge", alpha=0.5) +
    stat_pointinterval(.width = c(0.50, 0.9),position="dodge") +
    labs(x = "", y = "Perceived liking", title = "") + 
    scale_fill_manual(values=c("-1" = pal[1],"1"=pal[3]), guide="none") +
    scale_color_manual( values=c("-1" = pal[1],"1"=pal[3]), guide="none") +
    scale_x_discrete(labels=c("Tasteless", "Milkshake")) +
    scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(0,100, by = 20)), limits = c(-0.5,100.5)) +
    theme_bw()


plt_bayes_html = pp + html_theme 
plt_bayes = pp + averaged_theme 


figureHED <- ggarrange(plt_bayes_html, HDI_HED,
                    labels = c("A", "B"),
                    ncol = 2)
#figureHED


# AVERAGED EFFECT
dfH <- summarySEwithin(HED.means,
                       measurevar = "perceived_liking",
                       withinvars = c("condition","session"),
                       idvar = "id")

dfH$cond <- ifelse(dfH$condition == "1", -0.25, 0.25)
HED.means$cond <- ifelse(HED.means$condition == "1", -0.25, 0.25)
HED.means <- HED.means %>% mutate(condjit = jitter(as.numeric(cond), 0.3),
                                  grouping = interaction(id, cond))


pp <- ggplot(HED.means, aes(x = cond, y = perceived_liking,
                            fill = as.factor(condition), color = as.factor(condition))) +
  geom_point(data = dfH, alpha = 0.5) +
  geom_line(aes(x = condjit, group = id, y = perceived_liking), alpha = .3, size = 0.5, color = 'gray') +
  geom_flat_violin(scale = "count", trim = FALSE, alpha = .2, aes(fill = as.factor(condition), color = NA))+
  geom_point(aes(x = condjit, shape = as.factor(intervention)), alpha = .3,) +
  geom_crossbar(data = dfH, aes(y = perceived_liking, ymin=perceived_liking-se, ymax=perceived_liking+se), width = 0.2 , alpha = 0.1)+
  ylab('Perceived liking') +
  xlab('Taste') +
  scale_y_continuous(expand = c(0, 0), breaks = c(seq.int(0,100, by = 20)), limits = c(-0.5,100.5)) +
  scale_x_continuous(labels=c("Pleasant", "Neutral"),breaks = c(-.25,.25), limits = c(-.5,.5)) +
  scale_fill_manual(values=c("1"= pal[3], "-1"=pal[1]), guide = 'none') +
  scale_color_manual(values=c("1"=pal[3], "-1"=pal[1]), guide = 'none') +
  scale_shape_manual(name="Intervention", labels=c("Placebo", "Liraglutide"), values = c(1, 2)) +
  theme_bw()+ facet_wrap(~session, labeller=labeller(session = labels))


pp + html_theme +theme(legend.position=c(.96,.94))
A) Posterior distribution by group. B) Highest density interval (90% HDI) for the interaction condition by session

  1. Posterior distribution by group. B) Highest density interval (90% HDI) for the interaction condition by session 

plt = pp + averaged_theme +theme(legend.position=c(.96,.94))
# PLOT OVERTIME


pp <- ggplot(HED.p, aes(x = as.numeric(trial), y = perceived_liking,
                        color = condition,
                        fill  = condition))+
    geom_point(data = HED.group, aes(shape = intervention, color = condition), alpha = 0.3) +
  geom_line(alpha = .5, size = 1, show.legend = F) +
  geom_ribbon(aes(ymax = perceived_liking + se, ymin = perceived_liking - se),  alpha=0.4) +
  geom_point() +
  ylab('Perceived Liking')+
  xlab('Trial')+
  scale_color_manual(labels = c('Pleasant', 'Neutral'), name = "",
                     values = c( "1" =pal[3], '-1' =pal[1])) +
  scale_fill_manual(labels = c('Pleasant', 'Neutral'), name = "",
                    values = c( "1" =pal[3], '-1'=pal[1])) +
  scale_y_continuous(expand = c(0, 0),  limits = c(0,100),  breaks=c(seq.int(0,100, by = 20))) +
    #scale_x_continuous(expand = c(0, 0),  limits = c(0,21),  breaks=c(seq.int(1,21, by = 2))) +
  guides(color=guide_legend(override.aes=list(fill=c(pal[3], pal[1]), color=c(pal[3], pal[1]))))+
    scale_shape_manual(name="Group", labels=c("Placebo", "Liraglutide"), values = c(1, 2, 18)) +
  theme_bw() +
  facet_wrap(~session, labeller=labeller(session =labels))


pp + html_theme + theme(strip.background = element_rect(fill="white"), legend.justification = "top", legend.position="right", axis.text.x = element_text(size = 14))
Perceived liking for each condition and group over trials.

Perceived liking for each condition and group over trials. 

plt = pp + averaged_theme + theme(strip.background = element_rect(fill="white"), axis.text.x = element_text(size = 14), legend.justification = "top", legend.position="right")



##We see here the base difference between Placebo VS Liraglutide BUT it's the same for both sessions so .. no effect of the treatment really XXXX

# bmod_time = update(bmod_full,  ~ .+as.factor(trialxcondition)) #evaluate interaction
# 
# df =  bmod_time %>%
#      emmeans(~ condition:session:trialxcondition) 
# 
# 
# 
# bayes_plt = as_tibble(df) %>%
#     ggplot(aes(x = trialxcondition, y = emmean,  fill =as.factor(condition))) +
#     #geom_point(position = "jitter") +
#     geom_smooth(method=loess, se = T,inherit.aes = T) + 
#     #stat_slab(.width = c(0.50, 0.95),position="dodge", alpha=0.5) +
#     #stat_pointinterval(.width = c(0.50, 0.95),position="dodge") +
#     ylab('Mobilized effort (a.u.) \n Baseline corrected')+
#     xlab('')+
#     scale_fill_manual(values=c("1" = pal[2],"-1"=pal[1]), guide="none") +
#     scale_color_manual( values=c("1" = pal[2],"-1"=pal[1]), guide="none") +
#     #scale_x_discrete(labels=c("CS+", "CS-")) +
#     theme_bw() + facet_wrap(~session, labeller=labeller(session =labels))



cairo_pdf('figures/Figure_HED_Lira.pdf')
print(plt)
dev.off()

cairo_pdf('figures/Figure_HED_Lira_Bayes.pdf')
print(figureHED)
dev.off()

This part are not run here for different reasons (time it take and need to be on the cluster)

cwd=pwd
cd /home/OBIWAN/LIRA/fMRI/HED_GLM

Packages

Thanks to all the package’s authors and developers!

report::report_packages()
##   - repro (version 0.1.0; Aaron Peikert, Andreas Brandmaier and Caspar van Lissa, 2021)
##   - ggpubr (version 0.4.0; Alboukadel Kassambara, 2020)
##   - papaja (version 0.1.0.9997; Aust, Barth, 2020)
##   - tinylabels (version 0.2.1; Barth, 2021)
##   - cowplot (version 1.1.1; Claus Wilke, 2020)
##   - apaTables (version 2.0.8; David Stanley, 2021)
##   - Rcpp (version 1.0.6; Dirk Eddelbuettel and Romain Francois, 2011)
##   - Matrix (version 1.3.2; Douglas Bates and Martin Maechler, 2021)
##   - lme4 (version 1.1.27.1; Douglas Bates et al., 2015)
##   - XML (version 3.99.0.6; Duncan Temple Lang, 2021)
##   - intmed (version 0.1.2; Gary Chan, 2020)
##   - pander (version 0.6.4; Gergely Daróczi and Roman Tsegelskyi, 2021)
##   - ggplot2 (version 3.3.4; Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.)
##   - plyr (version 1.8.6; Hadley Wickham, 2011)
##   - stringr (version 1.4.0; Hadley Wickham, 2019)
##   - tidyr (version 1.1.3; Hadley Wickham, 2021)
##   - scales (version 1.1.1; Hadley Wickham and Dana Seidel, 2020)
##   - usethis (version 2.0.1; Hadley Wickham and Jennifer Bryan, 2021)
##   - devtools (version 2.4.2; Hadley Wickham, Jim Hester and Winston Chang, 2021)
##   - dplyr (version 1.0.7; Hadley Wickham et al., 2021)
##   - kableExtra (version 1.3.4; Hao Zhu, 2021)
##   - afex (version 0.28.1; Henrik Singmann et al., 2021)
##   - ggthemes (version 4.2.4; Jeffrey Arnold, 2021)
##   - optimx (version 2021.6.12; John Nash, Ravi Varadhan, 2011)
##   - cmdstanr (version 0.4.0.9000; Jonah Gabry and Rok Češnovar, 2021)
##   - JWileymisc (version 1.2.0; Joshua Wiley, 2020)
##   - tidybayes (version 2.3.1; Kay M, 2020)
##   - MBESS (version 4.8.0; Ken Kelley, 2020)
##   - tibble (version 3.1.2; Kirill Müller and Hadley Wickham, 2021)
##   - rlist (version 0.4.6.1; Kun Ren, 2016)
##   - sjPlot (version 2.8.8; Lüdecke D, 2021)
##   - bayestestR (version 0.10.0; Makowski et al., 2019)
##   - coda (version 0.19.4; Martyn Plummer et al., 2006)
##   - caret (version 6.0.88; Max Kuhn, 2021)
##   - lspline (version 1.0.0; Michal Bojanowski, 2017)
##   - brms (version 2.15.0; Paul-Christian Bürkner, 2017)
##   - R (version 4.0.4; R Core Team, 2021)
##   - psych (version 2.1.6; Revelle, 2021)
##   - BayesFactor (version 0.9.12.4.2; Richard Morey and Jeffrey Rouder, 2018)
##   - emmeans (version 1.6.1; Russell Lenth, 2021)
##   - Rmisc (version 1.5; Ryan Hope, 2013)
##   - janitor (version 2.1.0; Sam Firke, 2021)
##   - lattice (version 0.20.41; Sarkar, Deepayan, 2008)
##   - corrplot (version 0.89; Taiyun Wei and Viliam Simko, 2021)
##   - pbkrtest (version 0.5.1; Ulrich Halekoh, Søren Højsgaard, 2014)
##   - knitr (version 1.33; Yihui Xie, 2021)
# aaronpeikert/repro@devel #repro, crsh/papaja@devel,knitr, afex #aov_car/lmer, emmeans #emmeans, parallel #detectCores, tidyverse ##ggplot/dplyr/plyr/tidyr, car #densityPlot, lspline # lspline, JWileymisc #egtable, kableExtra #kable_styling, glmnet #cv.glmnet, ggthemes #theme_fivethirtyeight,  MBESS #pes_ci, sjPlot #tab_model/plot_model, Rmisc #SummarySEwithin, janitor #row_to_names, rlist #list.clean, sessioninfo #sessioninfo, stringr #str_list, XML #readHTMLtable, bayestestR #bayesfactor_models, cmdstanr #multithead, ggpubr::ggarrange #   Arrange Multiple ggplots, tidybayes #spread draws]